Electrostatic fluctuations in cavities within polar liquids and thermodynamics of polar 

solvation 
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We present the results of numerical simulations of fluctuations of the electrostatic potential and 
electric field inside cavities created in the fluid of dipolar hard spheres. We found that the ther- 
modynamics of polar solvation dramatically changes its regime when the cavity size becomes about 
4-5 times larger than the size of the liquid particle. The range of small cavities can be reasonably 
understood within the framework of current solvation models. On the contrary, the regime of large 
cavities is characterized by a significant softening of the cavity interface resulting in a decay of the 
fluctuation variances with the cavity size much faster than anticipated by both the continuum elec- 
trostatics and microscopic theories. For instance, the variance of potential decays with the cavity 
size Ro approximately as 1/Rq~^ instead of the 1/Ro scaling expected from standard electrostat- 
ics. Our results suggest that cores of non-polar molecular assemblies in polar liquids lose solvation 
strength much faster than is traditionally anticipated. 
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I. INTRODUCTION 

Solvation represents the change in the free energy when 
a usually molecular object is inserted into a condensed- 
phase environment. Since a significant part of chemistry 
and all life processes happen in liquid solutions, the tradi- 
tional focus has been on solvation in liquids, polar liquids 
in particular. The heterogeneous problem of solvation is 
probably as complex as the theory of liquids itself and is 
hunted by the same basic issues making the quantitative 
description of liquids so hard. There are two dominating 
and mutually compensating contributions to the free en- 
ergy of solvation: the positive free energy of creating a 
cavity (empty space) for a molecule to be inserted and a 
negative stabilization energy from short range (van der 
Waals) and long-range (electrostatic) forces The pos- 
itive cavity free energy is normally significantly compen- 
sated by the negative stabilization free energy resulting 
in the overall solvation free energy, a situation akin to 
the competition between repulsive and attractive forces 
in equilibrium liquids 0]- 

The present study is devoted to electrostatic solva- 
tion, i.e. the free energy arising from the electrostatic 
interactions between the charge distribution of the so- 
lute with the charge distribution of the liquid solvent. 
The charge distribution within molecular solutes is of- 
ten modeled by atomic partial charges efficiently used in 
force fields of numerical simulations. On the contrary, the 
charge distribution of the solvent molecules is often well 
represented by molecular multipoles following the well- 
established tradition of classical electrostatics 0] and di- 
electric theory Q . Extensions to models utilizing atomic 
charges are also possible as used in numerical simulations 
Q and interaction-site models of molecular liquids [2^]. 

Electrostatic solvation is believed to be well- 
understood. Following Born fd] and Onsager [7], the 
problem is traditionally recast in terms of continuum 
electrostatics where the electrostatic free energy is sought 
for the solute charges inserted into a dielectric cavity. 



This approach has been extensively tested against the ex- 
perimental database of solvation of small ions and neutral 
molecules in polar molecular liquids [8|]. Despite some 
inconsistencies, the formalism can be easily incorporated 
into quantum calculations and can even be quantitative 
once the dielectric cavity is properly parametrized. 

There are however still some fundamental issues that 
cannot be addressed within elecrostatic models. The so- 
lution of the Poisson equation in dielectric media is es- 
sentially a boundary condition problem in which the as- 
sumptions tacitly made by the material Maxwell's equa- 
tions about the structure of the dielectric interface are 
essential for the solution. The standard electrostatics 
assumes abrupt discontinuity of the dipolar polarization 
at the dielectric surface. This boundary condition cre- 
ates surface charge ^ which is ultimately responsible 
for the electrostatic potential within the dielectric cavity. 
Whether interfaces of real polar liquids match the as- 
sumption of abrupt discontinuity of the bulk polarization 
is an open question. For instance, the electric field within 
a cavity in a polar liquid was found to be much different 
from the prediction of standard electrostatics up to the 
cavity size of a mesoscale dimension [10]. 

A new additional piece of evidence comes from stud- 
ies of hydrophobic solvation essential for colloid sta- 
bility, biopolymcr folding, and formation of biological 
supramolecular structures [ll|j US] • It was found that sol- 
vation of non-polar solutes changes dramatically in char- 
acter at the length of about 1 nm, which is about three 
molecular diameters for aqueous solvation [T3| . Solvation 
of solutes larger than this characteristic length was found 
to be dominated by surface effects, i.e. the structure 
of water at the hydrophobic interface. Weak dewetting 
i.e. a substantial decrease of the water density at 
the interface compared to the bulk water, was found to be 
a central part of solvation of large hydrophobic solutes. 

Given the current interest in solvation at mesoscale 
0' [i3K to a large extent driven by biological applica- 
tions [jj, we address here the problem of electrostatic 
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solvation of solutes significantly larger than have been 
mostly studied so far. Our study is driven by the question 
whether the change in the solvation character established 
for hydrophobic solutes is reflected in an equally dra- 
matic change in the character of electrostatic solvation. 
The fact that the properties of a polar liquid interface 
are inconsistent with the assumptions of Maxwell's elec- 
trostatics [10] points to the possibility of a new solution 
once the size of the solute exceeds some critical dimen- 
sion. This is indeed the result we report here. 

We have found from numerical simulations that the 
scaling of the fluctuations of the electrostatic poten- 
tial and electric field with the cavity radius is consis- 
tent with the expectations of electrostatics (qualitatively) 
and molecular solvation models (quantitatively) for small 
solutes, but changes dramatically at approximately the 
same solute/solvent size ratio as observed for hydropho- 
bic solvation. It turns out that the core of the solute be- 
comes non-polar with its growing size much faster than 
is normally anticipated. We will start with formulating 
the general results of the Gaussian solvation thermody- 
namics and discuss the outcome of computer simulations 
next. 



II. THERMODYNAMICS OF ELECTROSTATIC 
SOLVATION 

By definition, the chemical potential of electrostatic 
solvation is given by the ratio of two partition functions: 
the one which includes the electrostatic solute-solvent 
potential Vqs and the one which is based on the non- 
electrostatic solute-solvent interactions and the interac- 
tions between the solvent particles. All these latter in- 
teractions are incorporated in the Hamiltonian Hq. The 
relation for /ios is then 



-/3mo=(/3) 



where 



Q(/3) 



dr. 



(1) 



(2) 



Here, we use the subscript "0" for the solute and the 
subscript "s" for the solvent, dT denotes integration over 
the system phase space, and (3 is the inverse temperature. 
Equation ([T]) can be conveniently re-written in terms of 
the product of the Boltzmann distribution of finding the 
solute-solvent energy e = Vbs and the probability density 



where 



P(e,/3) = Q(/3)'i / d{e-Vos)e 



-/3-f/o 



dr. 



(3) 



(4) 



Equation ([3]) is exact and it states that all the thermo- 
dynamic information required to understand electrostatic 
solvation is contained in the distribution of fluctuations 
of the interaction energy e = Vqs produced by the solvent 
which is actually not polarized by this potential; Vqs = 
for the Hamiltonian Hq. 

The approximation that we will adopt in our formal- 
ism, which is supported by our present simulations and 
data from other groups [Hi [13, [2l| , is to assume that the 
distribution function P(e, /3) is a Gaussian function with 
zero average 



P(e, /?) oc exp 



(5) 



The approximation of zero average is the reflection of the 
fact that no specific orientation of the solvent dipoles is 
created around a non-polar solute. This app roximation 
is not necessarily always correct [22l . l23j |. but is in- 
significant for most of our development since a non-zero 
average, if it exists, can always be incorporated in a lin- 
ear shift of e. What is the most significant property for 
our analysis is the magnitude and the temperature de- 
pendence of the Gaussian width cr^ (/3) . 

Within the Gaussian approximation for the electro- 
static fluctuations around a non-polar solute the ther- 
modynamics of solvation gains a simple and physically 
transparent form. The chemical potential of solvation is 

fios = ~{/3/2)a'{(3). (6) 

In addition, one can determine the energy e and entropy 
s of electrostatic solvation 



Ts = 



(7) 



In this equation, (Vbs) is the average solute-solvent elec- 
trostatic interaction energy when full solute-solvent in- 
teraction is turned on. From Eqs. ^ and ([7]), 

{Vos) - ~f3a'{f3). (8) 

The term Acgs in Eq. ([7]) determines the change in the 
interaction energy between the solvent molecules induced 
by electrostatic solute-solvent interaction. This energy 
term is identically equal to the corresponding contribu- 
tion to the solvation entropy, TAsss = Aess, so that Acss 
cancels out in the solvation chemical potential which is 
determined by solute-solvent interaction thermodynam- 
ics only [23, [2^. The term Acss can be calculated by 
either taking the derivative of the Gaussian width cr^(/3) 
or from a third-order correlation function 

Ae.. = = [(5V2){5Vl5H,),. (9) 

In Eq. ([H]), the average (. . . )o is over the ensemble of 
the non-polar solute in equilibrium with the solvent, col- 
lectively described by the Hamiltonian i?o- In addition, 
<5Vbs = Vbs - (Vbs)o and 5Hq = Hq - {Ho)o are devi- 
ations from the average values determined on the same 
unpolarized ensemble. 
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FIG. 1: Dielectric constant ta of the liquid of dipolar hard 
spheres vs the dipolar parameter {m*)'^ = f3m^/a"^; p* = 
0.8. The solid line represents the Fade approximation of the 
simulation data: €s{x) = (1 + aix + a2x'^)/{l + bix + b2x'^) 
with ai = 2.506, 02 = 3.057, fei = -0.180, 62 = -0.00865 and 
x = (m*)^. The dielectric constants were calculated from 
NVT MC simulations of the ho mog eneous liquid of dipolar 
hard spheres using the Neumann [20] correction for the cutoff 
of dipolar interactions treated by the reaction-field formalism. 



III. SIMULATIONS AND DATA ANALYSIS 

While the equations presented in Sec. |TT] are gener- 
ally applicable to an arbitrary solute, we will use nu- 
merical Monte Carlo (MC) simulations Q to determine 
the statistics of fluctuations produced in spherical cav- 
ities carved from a liquid of dipolar hard spheres (see 
Appendix for the description of the simulation protocol) . 
The fluid of dipolar hard spheres leaves out many im- 
portant properties of real liquids, most notably van der 
Waals forces and higher order multipoles. However, it 
allows a significant simplification of the solvation ther- 
modynamics since all physical properties of the solvent 
are expressed in terms of only two parameters, the re- 
duced density p* = pcr^ and the reduced dipole moment 
(m*)^ = 13m? / , where m is the dipole moment and 
cr is the diameter of the dipolar particles. Since the re- 
duced density is fixed to p* — 0.8 in our simulations, 
our results are fully defined in terms of two parameters: 
the reduced cavity radius Rq/u and the polarity param- 
eter (to*)^. The representation in terms of the dielectric 
constant can be easily achieved as well since these are 
well tabulated from our simulations as is shown in Fig.[T] 
The dielectric constants were calculated from Neumann's 
formalism [2^ as described in detail in Ref. 2lL 

We will also limit our consideration to two types of 
electrostatic multipoles most commonly studied in the- 
ories and applications of solvation, point ion and point 
dipole 0,0, [2l[. In both cases, the corresponding mul- 
tipole is placed at the center of the spherical cavity. 
The solute-solvent interaction potential is then given as 
Vbs — qo<j>s in the case of the ion and Vbs — —mo • Eg for 
the dipole. In these relations, and mg are the charge 
and dipole moment of the probe multipole and 0s and 
are, respectively, the potential and electric field produced 
by the solvent at the multipole position. 
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FIG. 2: Ai (a) and (b) vs the cavity radius Rq for 
(m*)^ = 0.5 (circles), 1.0 (squares), 2.0 (diamonds), and 3.0 
(up-triangles). The dashed line in (a) gives the result of Eq. 
(I17p for m* = 1.0. The dash-dotted and dashed lines in (b) 
shows the application of Eq. (|18|l at (m*)^ = 0.5 and 1.0, 
respectively. 

The main parameter entering the Gaussian model of 
solvation that we want to monitor is the Gaussian width 
CT^(/3). Since we want to deal with dimensionless quan- 
tities, we will in fact calculate the temperature reduced 
parameter 

T^p^<T\p)=P^{{5V^,fW (10) 

Since this parameter depends on the multipolar character 
of the solute, it is convenient to take this information out 
and consider the parameter A such that the temperature- 
reduced electrostatic energy of the solute is taken out as 
a multiplier 

r = u;A. (11) 

Here, the electric field of the multipole (charge or dipole) 
Eq is used to define the electrostatic energy 

w = {13/Stt) [ Eo{rfdr, (12) 
Jo. 

where the integral is taken over the solvent volume out- 
side the spherical cavity. 

The parameter w is equal to (3ql/{2Ro) for an ion and 
/3mo/(3i?g) for a dipole, where Rq is the cavity radius. 
Therefore, one can calculate the parameter A according 
to the following relations in case on ion (subscript "i") 
or dipolar (subscript "d" ) solutes 

A, = 2(3Ro{i6cbs)')o, 
A, ^ (3Ro{ms f)o- 
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Here, r^s = Ro/cr + 0.5 is the reduced distance of the 
closest approach of the hquid molecules to the cavity and 
y = (47r/9)/3TO^p is the standard density of dipoles in 
the dipolar liquid [1], p is the liquid number density. In 

addition, lQ^J{ros,p*) is the three-particle perturbation 
integral which is a function of the liquid density and r^s 
and i?cff(fos, P*) is the effective radius of the cavity 



Rj{ros,p*) 



dr 



(19) 



FIG. 3: {Ro/a)*A^ vs the cavity radius Ro for (m*)^ = 0.5 
(circles), 1.0 (squares), 2.0 (diamonds), and 3.0 (up-triangles) . 



Similarly we will introduce the reduced parameter A^s for 
the components of the internal energy and entropy arising 
from the alteration of the solvent-solvent interactions. 



Al=p^R^{{6^sfSHo)o, 
Ai = {p'Rl/2){{6Esr5Ho) 



(14) 



A few analytical results from standard electrostatics 
3 can be used as benchmarks in calculating and A^. 
The continuum electrostatics of Born @ and Onsager [7| 
equations gives the response functions A^^^ depending 
only on the dielectric constant of the dipolar liquid: 



A,: = 2 1 - - 



and 



6- 



2e, 



1 



(15) 



(16) 



In addition, several microscopic relations have been de- 
rived based on different formulations of the liquid-state 
theory. A closed-form equation for ion solvation is pro- 
vided by the Ornstein-Zernike integral equations for the 
ion-dipole mixture solved in the mean-spherical approx- 
imation (MSA) HI: 



A..^^fl-1 
Ro + Al \ es 



(17) 



In this equation, A^ = 3ct^/(1 -1-4^) is the correlation 
length of longitudinal polarization fluctuations of a dipo- 
lar liquid and ^ is the MSA polarity parameter [2^ . 

An analogous MSA solution exists for the mixture of 
dipolar particles of different size [30] which gives the pa- 
rameter Ad- Truncated perturbation expansions [3l| are 
however known to work better in this case with the result 
113,131 



A. = 6(|L 



y 



l + K(y,ros)y<7^4f/Rl 



(18) 



off 



In this equation g^J!) (r) is the hard-sphere distribution 
function of the liquid particles as a function of the dis- 
tance r to the cavity center. All functions Ros{rQs, p*), 
(^Os , P* ) 1 and k(?/, tqs) are given as analytical func- 
tions of the corresponding parameters in Ref. 132 . 



IV. RESULTS 

Our simulations have produced an unexpected result. 
We found that the scalings of electrostatic fluctuations 
and the corresponding chemical potentials with the cav- 
ity size do not follow the predictions of both the con- 
tinuum electrostatics and microscopic solvation models 
in case of large cavities. The results are shown in Fig. 
m As is seen, the parameter A^ decays much faster than 
the expected I/Rq scaling for all cavities greater than the 
size of the solvent particle. The large cavity scaling does 
not follow any universal law, but instead depends on the 
polarity (parameter to*) of the liquid (Fig. [3|). For the 
liquid polarities studied here, the large-cavity scaling of 
Ai is approximately l/i?g~^. Fluctuations of the electric 
field at the cavity center, representing dipole solvation, 
do not deviate that dramatically from the traditional ex- 
pectations, but the parameter A^ still decays to zero in- 
stead of leveling off as suggested by Eqs. (fT5|) and (fTSjl . 
In fact. Ad follows Eq. (fT5|) quite well up to the cavity 
size about 4-5 times larger than the liquid particle, but 
then starts to drop following qualitatively the trend seen 
for the potential fluctuations. Continuum electrostatics 
[Eq. (fTS)) ] fails both qualitatively and quantitatively for 
electrostatic fluctuations of both the potential and the 
electric field. 

There is a slight dependence of the variances on the 
number of particles in the simulation box. The variances 
extrapolated to A'^ — > oo from simulations done at various 
system sizes are listed in Table U in the Appendix. This 
dependence does not affect any qualitative conclusions 
we make here. Since extrapolation to iV — > oo creates a 
scatter of points, the results presented in Fig. [2] refer to 
a given system size only. 

With the dramatic failure of some very basic expecta- 
tions regarding electrostatic fluctuations, as is shown in 
Fig. [2l one wonders if the Gaussian approximation for the 
distribution of the electrostatic interaction energies fails 
for large cavities. We have tested this question by look- 
ing at the non-gaussianity parameter for both potential 
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FIG. 4: Non-gaussianity parameter So [Eq. (I20|l ] vs ts (a) 
and i?o (b). Points represent probe ions (circles) and probe 
dipoles (squares) for Ro/a — 2.5 (a) and (m*)^ = 0.5 (b). 



and field fiuctuations: 



(20) 



This parameter was found to be around zero, as expected 
for the Gaussian noise, within about 5% of the simula- 
tion uncertainties (Fig. [4]). The Gaussian approximation 
therefore seems reliable for our parameters database. 

In order to gain more insight into the origin of our 
observations, we have calculated two local parameters 
related to the orientational and density structure of the 
liquid/cavity interface. Figure shows the second-rank 
orientational order parameter of the permanent dipoles 
in the first solvation shell at the cavity surface: 



(21) 



Here, P2{x) is the second Legendre polynomial, fj — 
Yj/rj is the unit vector in the direction of the liquid 
particle j, and is the unit vector along its dipole mo- 
ment. The orientational order parameter shown in Fig. 
[5^ is calculated by limiting the distance r to liquid par- 
ticles residing in the cavity's first solvation shell where 
it indicates the existence of a preferential orientational 
order. The first-rank orientational parameter, based on 
the first-order Legendre polynomial, is identically zero 
thus implying that there is no net dipolar polarization 
at the cavity surface. This result is distinct from the 
water surface where water's large quadrupole moment is 
responsible for asymmetry 

As the cavity gets larger the solvent dipoles find it 
more energetically favorable to orient parallel to the in- 
terface, as was also observed for 2D dipolar liquids [33 |. 




FIG. 5: (a): the orientational order parameter vs the cavity 
size for different polarities of the solvent, (m*)^ — 0.5 (cir- 
cles), 1.0 (squares), 2.0 (diamonds), and 3.0 (up-triangles) . 
(b): contact value of the radial distribution function at 
Ri = Ro + a /2 vs the cavity radius. Shown are the results for 
different number of particles in the simulation box A'^ = 256 
(circles), 500 (squares), 864 (diamonds), 1372 (up-triangles), 
2048 (down-triangles), 2916 (stars), 4000 (pluses). Extrapo- 
lation to A'^ ^ oo is shown by bold solid line. The dashed 
lines connect the points. The thin solid line gives the contact 
value of the distribution function in the hard-spheres mixture 
from Ref. Isl. 
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FIG. 6: al3{{S(l>sf) (a) and a^l3{{SEsf) (b) vs the cavity size 
for probe charge and dipole located the distance a/2 from the 
cavity surface. The points refer to (m*)^ = 1.0 (circles), 2.0 
(squares), and 3.0 (diamonds); A'' = 1372. 



for water at cavity surfaces [35| and liquid-vapor inter- 
faces [9|, and for interfaces of dipolar liquids [s^ from 
density- functional calculations. However, this preferen- 
tial orientational order starts to dissolve with a further 
increase of the cavity size, after gaining maximum for the 
cavity about five times larger than the solvent particle. 
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This decay is related to the onset of softening of the first 
solvation shell indicated by the contact value of the pair 
cavity-solvent distribution function shown in Fig. [Hd- 

The contact value of the pair distribution function first 
rises as expected for a hard-sphere impurity in densely 
packed hard spheres [11] (solid hne in Fig. [3d), but then 
starts to drop. This drop appears at approximately the 
same value i?o/o' — 2 — 2.5 as both the downward turn of 
the orientational order parameter and the onset of devia- 
tion of the electric field fluctuations from the traditional 
predictions (Fig. [2]). We therefore can conclude that the 
observed change in the character of the electrostatic fluc- 
tuations is related to softening of the liquid/cavity inter- 
face, which also loosens the energetic push for a specific 
dipolar order. We note, however, that the peak of the 
distribution function stays at the closest-approach value 
Ri — Rq + a/2 and thus no dewetting of the cavity 
interface occurs. 

That the decay of the solvation energies is related to 
the softening of the interface is also seen from probing 
the fluctuations of the potential and field close to the 
cavity interface. Figure [6] shows the corresponding quan- 
tities for a point within the cavity kept one solvent ra- 
dius <t/2 away from the interface once the cavity size 
is increased. Again, simple electrostatic arguments sug- 
gest that the solvation energetics should approach that 
for a probe charge or dipole next to an infinite dielec- 
tric wall. Depending on how the dielectric interface is 
defined, by the cavity boundary or by the distance of 
the closest approach, continuum electrostatics predicts 
3 for af3{{64>s)^) the value between (cs - l)/{es + 1) and 
0.5(es — l)/(es -f f ). The observed dependence does seem 
to inflect into a plateau at the level consistent with this 
prediction at intermediate cavity size, but then starts to 
decay. This decay is however much more gentle than in 
Fig. [2] indicating that the area next to the interface is 
effectively stronger solvating than the part of the hollow 
space closer to the cavity center. 

The Gaussian approximation is a central part of our 
thermodynamic arguments and so we have done an ad- 
ditional test of its consistency also offering some deeper 
insights into the nature of electrostatic response func- 
tions. Since the chemical potential of solvation is given 
by the variance of the solute-solvent interaction potential 
[Eq. ©J, it becomes quadratic in a test multipole used to 
probe the electrostatic fluctuations. This result, known 
as the linear response approximation (2l| . suggests that 
the response function, obtained as the second derivative 
of in the corresponding multipole, does not depend 
any more on the magnitude of that multipole. ft also 
implies that Ai^d can be obtained from simulations of 
empty cavities but also from simulations involving ac- 
tual multipoles inside the cavity. The chemical potential 
of solvation and corresponding parameters ^ are then 
calculated from the average solute-solvent interaction en- 
ergy using Eq. ([8]). Since such simulations involving the 
probe charge are not straightforward due to the break- 
down of the system neutrality and the related difficulty of 
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FIG. 7: Response function Ad (a), contact value of the cavity- 
solvent pair distribution function gos{Rl) (b), and the orien- 
tational order parameter of the dipoles in the first solvation 
shell P2{Ri) (c) vs the magnitude of the solute dipole at the 
cavity center mo. Points are the results of MC simulations 
with Ro/a = 9.0, [m'f = 1.0, and = 2048. The lines 
connect the simulation points. 



using the Ewald sums [33, [Sg , we have done simulations 
of point dipoles of varying magnitude placed at the cav- 
ity's center. The results are shown in Fig.[7|for the cavity 
size above the threshold seen in Fig. [51 Rq/ct = 9.0, and 
the simulation box containing N = 2048 solvent particles. 
There we show the parameter calculated from (Vqs) at 
the varying magnitude of the solute dipole mg . Figures 
and [7J; also present the corresponding contact values 
of the cavity-solvent pair distribution function gos{Ri) 
and the orientational order parameter p2{Ri)- The re- 
sponse function A^; stays constant almost in the entire 
range of mg studied, starting to rise when the dipole in- 
side the cavity exceeds the solvent dipole by three orders 
of magnitude. This rise is a reflection of the change in the 
microscopic structure of the interface as the first solva- 
tion shell gets stiffer under the pull of the solute dipolar 
field and the first-shell dipoles start to reorient along the 
field of the solute dipole. The observed changes in the 
functions gQs{Ri) and p2{Ri) are, however, much greater 
than the corresponding change in A^ testifying to the col- 
lective nature of the solvent dipolar response effectively 
depressing changes in the microscopic structure of the 
first solvation shell. 

In Fig. [8] we show the same data as in Fig. [7l but ob- 
tained at a much smaller cavity size Ro/a = 1.5. Here, 
the change in the local structure with increasing the so- 
lute dipole is more pronounced and A^ starts to show 
a dependence on the magnitude of the probe dipole sig- 
nalling the appearance of nonlinear solvation effects. The 
variation in the response function is still mostly within 
10% and can be accounted for by nonlinear extensions of 



FIG. 8: Same as in Fig. [7] for Ro/a = 1.5. 



dipolar solvation models [27[. We need to stress, how- 
ever, that the Gaussian approximation appears to be ro- 
bust for large cavities which are of main interest for us 
here. 

We next turn to the dependence of the cavity response 
functions on the liquid polarity. Figure [9^ shows the 
dependence of on the solvent dipole moment. For 
a small cavity size, when the standard scaling with the 
cavity size is expected to apply, the dependence of 
on polarity does not show a saturation predicted by con- 
tinuum electrostatics [Eq. ([15])]. This saturation appears 
for a slightly larger cavity, but, as seen for a still larger 
cavity, it is simply en route to become a decreasing func- 
tion of polarity for the largest cavities studied here. We 
can therefore conclude that there is no range of param- 
eters where both the size scaling and the dependence on 
polarity predicted by the continuum electrostatics for the 
potential fluctuations are satisfied even at the qualitative 
level, not to mention the fact that the predicted values 
are significantly off. 

The saturation predicted by the Onsager equation for 
dipole solvation [Eq. is never reached. In contrast 

to the potential fiuctuations, the variance of the field is 
a uniformly increasing function with increasing solvent 
dipole for all cavity sizes studied here. A similar trend, 
for a narrower range of parameters, was previously ob- 
served by us [395 and it manifests itself in the solvation 
dynamics uniformly slower than continuum predictions 
[4q |. The results for A^ from Eq. arc shown by the 
solid lines in Fig.[5]3. As expected, there is a good agree- 
ment with the simulations for small cavities, but then the 
theory fails when the regime of solvation changes and A^; 
turns downward as in Fig. [51 

In Fig. [TU] we show the results of calculations of the 
solvent-solvent component of the solvation entropy [Eq. 
([7])]. Within the Gaussian approximation, the ratio of 



FIG. 9: Ai (a) and A^j (b) as functions of (m*)^ for 
i?o/(T = 0.5 (circles), 1.0 (squares), 1.5 (diamonds), and 6.0 
(up-triangles) . The solid line in (a) shows the result of using 
Eq. p7|l &t Ro/a = 1.0. The solid lines in (b) show the result 
of Eq. ([HJ for Ro/a = 0.5, 1.0, and 6.0 (from down up); the 
dashed lines in (a) and (b) connect the points. The data for 
Ai at iio/o" = 6.0 (up-triangles) in (a) have been multiplied 
by a factor of 5 to bring them to the scale of the plot. The 
simulation points were obtained at = 1372 dipolar hard 
spheres in the box. 



the solvent-solvent component of the solvation entropy, 
Tsss — Acss, and the solute-solvent component, Tsqs — 
— /3tT^(/3)/2, is given as the ratio of the corresponding 
reduced response functions 



i,d 



2A^ 



(22) 



As is seen, for both the ionic and dipole solvation, there 
is a compensation between the ordering of the solvent by 
the solute, expressed by always negative sqs, and the dis- 
ordering of the solvent structure, expressed by positive 
Sss- This compensation is however far from complete, 
in contrast to a much stronger compensation found for 
aqueous solvation [4l[. The overall entropy of electro- 
static solvation is therefore negative. Since the parameter 
Xs in Eq. ([^ depends weakly on the cavity size, the dra- 
matic change in the character of solvation found here for 
(T^(/3) will be reflected in both the enthalpy and entropy 
of electrostatic solvation which are often more accessible 
experimentally than solvation free energies. Very httle 
is currently known about the magnitude of Xs |42l |. in 
particular for large solutes. Our recent MD simulations 
of the redox entropy of metalloprotein plastocyanin (43j 
have produced Xs — 0.4 {Rq/ct ~ 5.8), although it is not 
clear if the Gaussian approximation is applicable to the 
protein electrostatics. 
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FIG. 10: The ratio of the solute-solvent and solvent-solvent 
components of the solvation entropy vs the cavity radius cal- 
culated for charge and dipole probe multipoles. 



V. DISCUSSION 

In this paper we have suggested to study polar solva- 
tion by using Eq. ([3]) which states that all the informa- 
tion required to calculate the solvation thermodynamics 
is contained in the distribution of electrostatic interaction 
energies around a fictitious solute with the solute-solvent 
electrostatic coupling switched off. This equation is exact 
and the approximation adopted here is that the distribu- 
tion function P(e) can be approximated by a Gaussian. 
The distribution P(e) can generally be written as 

P(e) oc exp [/3w(e)] (23) 

and then the integral over e in Eq. ([3]) can be taken by the 
steepest descent around the stationary point eo defined 
by the condition w'(eo) = 1. The Gaussian approxima- 
tion is then equivalent to assuming all the terms except 
the linear one can be dropped from the series expansion 
of w'(e) in powers of (e — eq). 

Our simulations have not identified any significant 
deviations from non-gaussianity. Extensive simulations 
done with ionic and dipolar solutes over the last decades 
0, [13, HH, have also resulted in the conclusion that 
the Gaussian picture is an accurate one implying that 
P(e) is globally a Gaussian function. However, one can 
argue that we could not sample sufficiently around eo 
and thus cannot assess the deviations from Gaussian- 
ity. While that might be true for strong solute-solvent 
interactions, for which a significant data-base pointing 
otherwise exists [sl], energy eo is expected to decrease 
with increasing the cavity size and the Gaussian approx- 
imation is expected to become increasingly accurate (as 
indeed seen from comparing Figs. [7] and [5]). However, it 
is in this range of large cavities, almost completely ne- 
glected in previous studies of electrostatic solvation, that 
we found the most dramatic deviations from the tradi- 
tional expectations. 

The main funding of this study is that electrostatic 
solvation by polar liquids changes its regime at the size 
of the cavity about 4-5 times larger than the size of the 
solvent particle. The regime of small cavities can be rea- 
sonably understood with molecular solvation models and 
in particular the results for the electric field fiuctuations 



(probe dipole) are in a very good quantitative agreement 
with the results of perturbation solvation models. The 
regime of large cavities is dramatically different and can- 
not be described by the models traditionally employed 
for solvation problems. 

What we have observed here is a dramatic decay of the 
solvation strength in the middle of the cavity, much faster 
than expected from both the continuum electrostatics 
and microscopic solvation models. For instance, the vari- 
ance of the electrostatic potential decays as 1/Rq~^ in- 
stead of the expected l/i?o scaling. The core of a growing 
hollow cavity thus becomes non-polar much faster than 
previously anticipated. What it practically means is that 
there is very little solvation stabilization for charges in- 
side a large mesoscale object. This might be a reason 
why natural systems requiring hydration of large molecu- 
lar assemblies (proteins, etc.) rely on solvation of surface 
charges for which much slower decay of solvating power 
due to softening of the interface was found here. In appli- 
cation to the problem of protein folding, this observation 
implies a very strong driving force for placing ionized 
residues and cofactors stabilizing protein solvation closer 
to the interface. 
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APPENDIX A: SIMULATION PROTOCOL AND 
RESULTS 

MC simulations were performed with the standard 
NVT Metropolis algorithm. The initial configuration was 
constructed starting from an fee lattice of liquid hard- 
spheres of diameter a and density p* = 0.8. The hard- 
sphere solute/cavity was then "grown" in the center of 
the simulation box by increasing the initial cavity di- 
ameter of 0.5(7 with 0.002ct increment, adjusting a to 
ensure constant density, and moving the solvent parti- 
cles according to the Metropolis algorithm. After the 
solute/cavity was constructed, the initial configuration 
was created from 10^ — 10^ parallel steps (using Open- 
MPI) producing different initial configurations for each 
processor. The subsequent runs were then carried out on 
each processor separately thus minimizing interprocessor 
communications. To guarantee the Markovian statistics, 
the random number generators used in the MC moves 
were seeded independently between the processors. This 
implementation has resulted in a linear scaling of the 
program output with the number of processors. The pro- 
duction runs of (1 — 5) X 10^ steps were performed on 10 
processors per (m*)^ per cavity size. 
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The simulation protocol employed the minimum im- 
age convention and the reaction-field correction |5|] for 
the cut-off of dipolar interactions at one-half of the cu- 
bic simulation box. Ewald sums were also tested 
and gave results identical within simulation uncertain- 
ties. The reaction-field correction was preferred due to 
better performance. The dependence on the simulation 
box size was carefully checked in particular since growing 



cavity required larger number of liquid particles to elim- 
inate finite-size effects. The number of particles N was 
varied in the range TV = 108, 256, 500, 864, 1372, 2048, 
2916, and 4000 depending on the cavity size. The rep- 
resentative results for and listed in Table U were 
obtained by averaging over several simulation runs with 
different box sizes and also by extrapolating the plots of 
corresponding values vs 1/A^ to the iV — > oo limit. 
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TABLE I: Values of Ai and for (m*)^ = 0.5,1.0,2.0,3.0, e = 3.63,8.51,29.9,93.7, respectively. Extrapolations (ext.) 
were done with N = 108, 256, 500, 864, 1372, 2048, 2916, 4000 data when available, linearly fitting Ai,d vs. 1/N and taking the 
intercept. The system sizes used for the extrapolations are given in the footnotes. 
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Ai 


Ad 
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A. 
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3.0 

A. 


Ad 
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1372 
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0.498 


1.031 
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1 


,185 


1.126 




ext. ^ 
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0.310 


2.004 


0.251 


2, 


,388 


0, 


,246 


2.556 


6.0 


1372 


0.075 


0.675 


0.061 


1.084 


0.051 


1 


,483 


0, 


,046 


1.670 




ext.' 


0.298 


1.449 


0.296 


2.061 


0.236 


2, 


,484 


0, 


,233 


2.553 


7.0 


1372 


0.043 


0.478 


0.033 


0.789 


0.024 


1 


,124 


0, 


,021 


1.263 




ext.' 


0.192 


1.233 


0.183 


1.834 


0.133 


2, 


,203 


0, 


,125 


2.373 


8.0 


1372 


0.026 


0.338 


0.018 


0.581 


0.012 


0, 


,855 


0, 


,010 


0.952 




ext.s 


0.128 


1.003 


0.114 


1.540 


0.054 


1, 
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9.0 


1372 
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0.803 
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0.036 
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